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The ideal invariants present in the formalism of magnetohydrodynamics (MHD), i.e. global quantities that are 
conserved in the absence of sources and dissipative effects, play an important role in various theoretical and nu- 
merical studies of MHD turbulence. The fluxes of these ideal invariants represent separate channels that transfer 
the information across different scales in a turbulent system. Once a statistically stationary state of turbulence 
is reached, the amount of any ideal invariant quantity introduced in the system by a forcing mechanism equals 
the amount of the same quantity removed by the dissipative effects from the system. For highly developed tur- 
bulence, these two mechanisms act predominantly at different scales that are largely separated. Since the ideal 
invariant quantities cascade between scales, a constant flux is generated with great implication on the state of 
the system. Numerically, controlling the ideal invariant fluxes levels for a turbulent MHD system is important 
for the analysis of fundamental MHD turbulence properties. We propose a forcing mechanism that controls the 
three ideal invariants of MHD turbulence: the total energy, the cross-helicity and the magnetic helicity. This 
forcing is implemented in the freely available TURBO solver, that is also briefly presented. 



L INTRODUCTION 

The motion of a fluid is described mathematically by a nonlinear evolution equation, for which no 
general analytical solution exists. Except for laminar flows in very simple geometries, even the use of 
a perturbative approach is impossible. This known problem, caused by the coupling of different scales 
in a flow by the nonlinear terms has triggered the use of numerical methods in the study of fluid flows. 
Usually, solutions of a flow are obtained using a numerical solver, in an incremental manner, starting 
from a set of initial conditions and a prescribed set of boundary conditions. For a typical engineering 
problem, the correct description of the physical properties and geometry of the boundary conditions is 
crucial, as they are responsible for an entire class of instabilities introduced in the flow which in turn 
characterizes the subsequent motion. However, for the study of fully developed turbulence, where a 
huge range of scales are present and a clear separation exists between the geometry dependent large 
scales and the universal, dissipative small scales, periodic boundary conditions represent the preferred 
choice. In this situation, the driving instabilities are introduced in the form of an external force that acts 
only at a particular scale, usually a large one. The use of periodic boundaries conditions for a flow has 
the advantage of enabling us to easily translate the motion problem in terms of Fourier modes which 
can be solved numerically using spectral solvers. For this purpose we make use of the TURBO code, 
which we will briefly introduce in Section II of this paper. 

When the fluid is electrically conducting like in the case of a plasma or of a liquid metal, the mo- 
mentum balance equation is influenced by the Lorentz force and the number of non-linear terms in- 
creases compared to the case of fluid turbulence. This coupling of the Navier-Stokes equation with the 
Maxwell reduced equations, which leads to the magnetohydrodynamic (MHD) equations, increases the 



* bogdan.teaca@epfl.ch, Tel: +41-21-693.43.05 
T clalescl@jhu.edu 

* bknaepen@ulb.ac.be 



§ dcarati@ulb.ac.be 



2 



complexity of the interaction between scales. Compared to hydrodynamical flows, the ideal invariants 
change in the case of MHD, as the interplay between the velocity and the magnetic fields needs to be 
considered. The importance of the MHD ideal invariants will be discussed in Section III. 

MHD turbulence, considered as an initial value problem with periodic boundary conditions, tends 
to decay in absence of driving mechanism. To maintain a stationary turbulent state, it is needed to add 
sources in the evolution equations that mimic the effects of the various instabilities that may appear 
in the large, geometry dependent scales of a realistic system. One way of achieving this would be 
through a constrain on the large scale value of one or both fields. This is often implemented in spectral 
space by freezing the value of a limited set of Fourier modes characterizing the largest scales in the 
system [10, 12]. Alternatively, an external artificial force can be used [1, 15]. We propose such a force 
in Section IV which controls the injection rate of the three ideal invariants of MHD turbulence: the 
total energy, the cross-helicity and the magnetic helicity. 

II. TURBO SOLVER 

The TURBO code 1 is designed to solve the equations for an incompressible fluid in a three di- 
mensional slab geometry with periodic boundary conditions in the three directions. The real space 
representation of the incompressible MHD equations under the influence of an external force and in 



the presence of a constant magnetic field Bo are written as, 

— = u • Vu + b • Vb + B • Vb + iA7 2 u + f" - Vp + f c , (1) 

at 

dh 

— = u ■ Vb + b • Vu + B • Vu + r?V 2 b + f 6 , (2) 
V-u = 0,V-b = 0, (3) 



where u = u(x, t) is the fluid velocity field, b = b(x,t) is the magnetic field expressed in Alfven 
units and p = p(x, t) is the total, hydrodynamic + magnetic, pressure field divided by the constant 
mass density. Due to the incompressibility condition, the pressure p is not an independent variable and 
can be formally eliminated by solving the Poisson equation, 

V 2 p = -Vu : Vu + Vb : Vb . (4) 

The kinematic fluid viscosity is v and the magnetic diffusivity is r\. The divergence free external force 
fields f" = f"(x,i) and f b = f b (x, t) act on the velocity and magnetic fields, respectively. The two 
forces are part of a forcing mechanism that imposes the injection rates of the MHD ideal invariant 
quantities. A kinetic only forcing method (f b = 0), used previously in the literature for similar studies 
[7, 15], can also be employed. In the velocity evolution equation, f" can be considered as divergence 
free since the pressure will enforce the incompressibility of the velocity field by eliminating any V • f " 
contribution of the force. On the other hand, f b must always be divergence free as a consistency 
condition for the magnetic field. If desired, a Coriolis force f c = l~i x u acting on the flow can also be 
considered. The force appears as result of a reference system rotation, with the angular velocity fi. In 
that case, the centrifugal acceleration that depends explicitly on the distance to the rotation axis can be 
lumped into the pressure term due to the incompressibility condition and the use of periodic boundary 
conditions is still appropriate. 

In addition to the MHD equation, the TURBO code can also solve the evolution equations for a set 
passive scalars fields: 

^ = _u • Vc Q + k q V 2 c q + a a ({cp}) (5) 



The code can be downloaded freely from: http:// 
aqua . ulb . ac . be/turbo 
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where c a = c a (x, t) are passive scalar(s), each of which is characterized by a diffusion coefficient K a . 
There is the possibility to include source or sink terms or even chemistry terms in the scalar equations 
through the function a a . 

In the TURBO code, the space discretization is based on a Fourier representation of the quantities 
of interest. For a given quantity, the physical Q and the spectral Q representations are related using the 
direct and the inverse discrete Fourier transforms 2 : 

Q(k) = ^£Q(x) e -* kx , (6) 

X 

Q(x)=^Q(k) e lkx . (7) 

k 

where N is the total number of modes in a given direction. In practice, the code allows the use of 
different numbers of modes in the three directions N x , N y and N z and different box sizes L x , L y 
and L z . Knowing the number of modes and the box size for each direction allows us to define the 
wavenumber space. Assuming a cubic box for simplicity of notations, the wavenumbers are defined 
as: 

k n = -j-n , (8) 

where n 6 [— N/2+1, N/2]. From the above definition, we see that the smallest nonzero wavenumber 
is fco = 2ir/L, which for the typical box length choice L = 2tt, becomes unity. The largest wavenum- 
ber accounted in a simulation (fc max = irN/ L) and the spacing between two neighbour grid points 
(A = L/N) are related as A = 7r/fc max . 

The main advantages of spectral methods are definitely the accuracy and the simplicity of the 
representation of the differentiation operator V, which reduces to a multiplication in spectral space, 
V u(x) — > iku(k). The weak point of this method is the mandatory use of periodic boundary condi- 
tions, a choice considered as hard coded in the TURBO solver. Spectral methods are thus not adequate 
for exploring very complex geometries but are extremely useful for investigating the fundamental prop- 
erties of turbulence. 

The time evolution is based on a modified Williamson, four-step, third-order low storage Runge- 
Kutta method [18]. Since the equations are discretized in space, it is desirable that the transfer of 
information be limited to neighbor grid point in one time step At. This implies the CLF criterium 
which simply states that the time step At has to be smaller than the time necessary for the wave with 
the largest propagation speed (determined on the velocity ||u|| max or on the Alfven velocity ||b|| max ) 
to propagate between the smallest distance present between two grid points (A). The linear terms 
are solved in an analytical manner for each mode k by performing an appropriate change of variables 
u — > u'(u, k, i/, f2) 3 and b — > b'(b, k, rf) before formally performing the Taylor series expansion. 
Although this will still give us a third-order accuracy in a global sense for the solution, the linear terms 
do not affect the value of the time step. 

Since the evolution equations are solved is spectral space, the value of the nonlinear terms for a mode 
k has to be computed. Determining the nonlinear terms directly in Fourier space, although possible, 
is prohibitive from a computational cost point of view, as they are represented by convolutions of 
modes. Instead, inverse Fourier transforms are used for the velocity and the magnetic field and the 
nonlinear terms are computed in real space before being Fourier transformed back to the spectral 
representation. In practice, the nonlinear terms appearing in the evolution equation for a velocity mode 
correspond to the divergence of a symmetric tensor UiUj — bibj, while the nonlinear term appearing in 
the evolution equation for a magnetic field mode appears as the divergence of an anti-symmetric tensor 



2 Numerically, a Fast Fourier Transform algorithm is em- 
ployed through the use of the FFTW libraries, [13], which 
can be downloaded at: http : / /www . f f tw . org/ 



3 The operator of this tensorial transformation reduces to a 
diagonal operator for CI = 0. 
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Uibj — biUj. This has an impact on the number of Fourier transforms that need to be perform: six for 
the symmetric tensor and only three for the anti-symmetric tensor, instead of nine transforms for each 
term. Furthermore, the trace of the tensor UiUj — fybj can be lumped into the pressure term, leading 
to another reduction of the computational efforts. Two orthogonal projections with respect to k of the 
resulting nonlinear terms ensure that the solenoidal conditions are enforced for both the velocity and 
magnetic fields. 

Computing the nonlinear terms in real space reduces the numerical method type from a purely spec- 
tral method to a pseudo-spectral one. This approach of computing the nonlinear terms, although much 
faster then computing the convolution, adds another problem known as aliasing. The aliasing error 
originates in the fact that a plane wave with a wavenumber k takes exactly the same values on the grid 
as a plane wave with a wavenumber k + Nko. When the nonlinearities are computed in real space, the 
aliasing error becomes a serious issue. For a one dimensional quantity q represented by modes with 
k — k x [(—N/2 + 1), N/2], its square q 2 has modes that correspond to k = fc x [(—N + 2),N]. 
As a consequence, the mode of q 2 that corresponds to k = ko(—N + 2) is undistinguishable from 
the mode k = 2ko. Two approaches, known as dealiasing methods, are considered to eliminate this 
difficulty [14]. 

The first method consists in assuming that only the modes of q with k = fco x [(— M/2 + 1), M/2] 
are non zero and imposing that the other modes remain zero after each nonlinear computation. The 
two-third method consists simply taking M = 2/3 N and keeping all the modes outside the range 
k = ko x [(—N/3 + 1), N/3] to zero. This method represents the simplest way of fully removing the 
aliasing error. 

The second method is based on the property that shifting the grid by a distance d results in a modifi- 
cation for the modes by a phase e lkd . For one-dimensional systems, the aliasing error can be removed 
exactly by computing the nonlinearity twice on two grids shifted respectively by A/2 and — A/2 and 
by summing the two computations. The contributions that do not lead to aliasing errors are unaffected 
by this procedure. In three dimensional systems, height evaluations with different shifts are needed 
to get an alias-free computation of the nonlinear terms. This is of course quite prohibitive and an 
approximation is made. The height shift computations of the nonlinear terms are done as part of the 
sub-step of a Runge-Kutta scheme. The phase-shift dealiasing method represents the reason behind 
the four step implementation for the third order Williamson Runge-Kutta time advancement method. 
This dealiasing method represents an approximative process of removing the aliasing errors. 

The TURBO code is parallelized through the use of MPI. The numerical cube is split in the y- 
direction in spectral space, which limits the maximum number of processors used for one run to N y . 
To improve the use of available parallel computing resources, an additional parallelization direction is 
made over the number of instances performed during each run, i.e. number of "cubes" solved at the 
same time that can exchange data among themselves. This last feature of the TURBO code allows the 
direct computation of ensemble average quantities. 

III. MHD IDEAL INVARIANTS 

In MHD turbulence, the three quadratic 4 ideal invariants have an important role in the dynamics of 
turbulence and the resulting cascades that appear for fully developed turbulence regimes [4, 9]. For 
ideal MHD, that is an inviscid flow with zero magnetic diffusivity, total energy, cross-helicity and 
magnetic-helicity are conserved in absence of forcing. Magnetic-helicity H m = (a • b) is a purely 
magnetic invariant, defined as the scalar product between the magnetic potential a and the magnetic 
field b = V x a, where (...) denotes volume average. Since it possesses an inverse cascade (transfer 
form small to the large scales [8]), magnetic-helicity plays an important role in dynamo phenomena, see 
[2, 3]. The other two quadratic invariants, total energy E = E U + E b and cross-helicity H c = (u • b) 



4 Other invariants may exist, but only quadratic invariants 
are robust enough to survive truncation due to numerical 
discretization [5]. 
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FIG. 1. Example of the real-space density of cross-helicity level for balanced and imbalanced MHD turbulence 
solved for 256 modes in each direction. For the balanced case p c w while for the imbalanced case p c « 0.6. 
The lower panels depict histograms of the cross-helicity level over the entire computational domain. 



should be studied together as they affect each-other. Since kinetic E u = (u • u) /2 and magnetic E b = 
(b • b)/2 energies are not conserved individually, the use of Elsasser variables z ± = u ± b might be 
more appropriate. In fact, in the Elsasser representation, the cross-helicity and total energy information 
is contained in the definition of two ideal invariants E + = (z + • z + )/4 and E~ = (z~ • z _ ) /4, known 
as pseudo-energy. The E ± ideal invariants are positively defined, which represents an advantage in 
spectral study. Since cross-helicity is related to the degree of alignement between the velocity and the 
magnetic field, the cross-helicity level defined as, 

P E ~ (E+ + E-) ' W 

provides global information regarding the alignment present in a system (p c G [— 1, 1]). Point-wise, the 
presence of a large cross-helicity level p c (x) = iJ c (x) /E(x) gives rise to the phenomena of nonlinear 
depletion in the evolution equations. Nonlinear depletion weakens the ability of the nonlinear terms to 
mix the flow compared to hydrodynamic case, which results in a different turbulent mixing time for 
MHD turbulence [19]. Moreover, it was shown in [11] that MHD turbulence is composed of zones 
of highly aligned (p c (x) ~ +1) and highly anti-aligned structures (p c (x) ~ —1) which are thought 
to develop naturally through the process of dynamical alignment [6]. Therefore, a non-zero value 
for the global parameter p c denotes a preference in the generation of one type of aligned structures. 
This situation is known as imbalanced MHD turbulence, Figure 1. It is interesting to note that the 
presence of cross-helicity does not necessarily affect the equipartition of energy between the u and b 
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fields. However, it does change the pseudo-energy levels of z + and z and therefore the strength of 
the co-propragating and contra-propagating Alfven waves that scatter on each-other. 

Although the ideal invariants are not conserved in the presence of viscosity and magnetic diffusivity, 
they are still redistributed between different scales without loss or gain through their nonlinear fluxes. 
For X, Y, Z standing in for the fields, the flux through a spherical surface in the Fourier space defined 
by a radius k c is defined as: 

n£ z (fc c ) = (Z(k) ■ vY(k|fc < k c ) ■ x(k|fc > k c )) . (io) 

The average represents here the sum over all the Fourier modes. For a turbulent state, the value of the 
fluxes in the inertial-inductive range are expected to be constant. In the Elsasser formalism, the values 
11+ _ and III , reached in the inertial-inductive range, enter in the phenomenological definition of 
pseudo-energy scaling laws [16]. 



IV. THE FORCING MECHANISM 



As it is known, a forcing mechanism can be used to reach different turbulent regimes. In this study, 
the injection rates of the nonlinear ideal invariants will be used as the control parameters. These 
injection rates can easily be computed from equations (1-2). For instance, the injection rates of total 
energy and of cross-helicity (H c = (9?{u(k) • b(— k)})) are given by 

= (f"(k) • u(-k)) + (f b (k) • b(-k)) = e u + e b = e , (11) 

/ 

= (f»(k) • b(-k)) + (f b (k) • u(-k)) = a u + a b = ea, (12) 

/ 

where, thanks to the Parseval theorem, the volume average (...) can be identified as the average over 
the number of modes. The parameters e u and e b represent the power injected by f " and f 6 respectively. 
Since the cross-helicity is bounded by the total energy, the sum of a u and <r b , which individually denote 
the cross-helicity injected by f " and f b respectively, has to fulfil the condition: — e < a u + a b < +e. 
As such, the cross-helicity parameter a is bounded in the interval [—1,1]. Selecting the force control 
parameters in such a way to fix e and a will enforce the dissipation level for the energy and cross- 
helicity, once the stationary regime is reached as shown on Figure 2. This behavior is true only for 
ideal invariant quantities. For example, selecting the injection level of kinetic helicity will not enforce 
the kinetic helicity dissipation for MHD turbulence as it would in a purely hydrodynamic flow. 

Numerically, we consider the forces f "(k) and f fc (k) to be local quantities in Fourier space, which 
act in the same manner on all the modes Nf within a wavenumber shell defined by the interval Sf — 
[fcj n f , fc sup ]. Since u, b and f u ' b are divergence free, we use a helical decomposition [17] for the 
definition of the force. The helical decomposition projects a vector a(k) on a complex basis h± : 

a(k) = a+(k)h + + a_(k)h_ , (13) 

where h± = ei x e2 ± iei and ei = (A X e2 ) / 1 1 A X ©2 1|« e2 = k/fc. The wave-vector A is taken to be 
arbitrary and non-parallel to k. We warn the reader that in this section, the ± lower indices refer to the 
helical basis h± , the vector projections on this basis and their contributions to the different scalar quan- 
tities of interest, while the upper indices denote, as usual, quantities in the Elsasser representation. The 
helical decomposition is very useful in ensuring zero divergence and in computing the curl operator. 
Indeed, the vectors h± are eigenmodes of the curl operator, ik x h± = ±fch±. As such the vorticity 
is now defined as uj± (k) = ±kii± (k) and the electric current has the form j± (k) = ±kb± (k). 
We choose the projection of the forces on h±, to have the form: 

f£(k) =a£(k)« ± (k)+/3£(k)6±(k) » (W) 
f£(k) =4(k)«±(k)+4(k)6±(k) , (15) 



dE 
~dt 
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FIG. 2. The total energy dissipation (D ) level (left) and and cross-helicity dissipation (D c ) (right) for two 
different a values. The plots are made for well resolved turbulence using a numerical resolution of 512 in each 
direction. 



if |k| 6 Sf and zero otherwise. Since the a's and /3's parameters that need to be determined are 
considered to be real, the forcing method presented here does not influence the phases of the fields, 
which ensures that no change is made in the type of turbulent structures present. Another way of 
injecting cross-helicity into the system would be achieved by imposing the alignment of u and b in the 
real space, which in turn would modify the phases of the fields and potentially the turbulence behavior. 

For a mode k, the kinetic energy (E u (k) = ^u(k) • u(— k)), magnetic energy {E b (k) = |b(k) • 
b(-k)), cross-helicity (H c {k) = 5R{u(k) • b(-k)}), kinetic-hehcity (H k (k) = 5R{u(k) • cb(-k)}) 

and the magnetic-helicity (H m (k) = 5R{b(k)-a(k)*} = —^di{h(k) • j(— k)})can be easily expressed 

k 

in terms of the helical decomposition. Knowing that for the helical decomposition, the energy injected 
per unit of time in the velocity equation and the magnetic equation are e u = e" + and e b = e b + + e b _ 
and that the cross-helicity injected per unit of time in the velocity equation and the magnetic equation 
are a u — c" + c" and a b = a b + + cr_, respectively, we can write the injection rates for the three ideal 
invariants and kinetic helicity per mode k in the forcing range as: 



dE(k) 
dt 
dH c (k) 

dt 
dH m {k) 
dt 

dH k (k) 
dt 



f N f 



1 [(el + el) + (e b + + e b _)] = ~e, 

TV/ 



/ N f 

- —-(e b -e b ) 

f N f 



(16) 
(17) 
(18) 
(19) 



To simplify the numerical implementation of the force, we make the following assumption: we 
consider a± = ae± and a b ± = ae b ± . In practice, these equalities impose that both forcing i u and 
f b are responsible for the same amount of cross helicity injection in the system. Moreover, they also 
impose that cross helicity is injected at the same rate in both the h + and h_ components of the velocity 
and magnetic fields. These assumptions, which could be reconsidered in the future, allow us to fix the 
eight real parameters a± (k), /?" (k), a b ± (k) and fi b ± (k) by giving only five control parameters, namely 
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the energy injection rates e±, e± and the cross-helicity parameter a, 



s _ e± eri?l(k) 2 - 2E b ± (k) 
a± (k) = iV^ i?|(k) 2 -4^(k)S|(k) ' (20) 



o-ff|(k) 2 - 2£|(k) 



^ (k) ^>iT^(k) 2 -4£^(k)^(k) ' (21) 



o-^(k) 2 - 2E b ± (k) 
N f Hl(k) 2 -4E|(k)E^(k) 



4 00 = ^r~^rr^:,, . (22) 



1 e5_ erffifk) 2 - 2K(k) 
tf b (LA — _± ±v ; ±v ; (?^n 

P± [ ' N f HI (k) 2 - 4£| (k)E b ± (k) ' K ' 

where the respective injection rates are assumed to be the same for all the Nf forced modes. Selecting 
a large number of forced modes ensures that no anisotropy effect is induced by the forcing mechanism. 
Because of the Cauchy-Schwarz inequality, we have the condition H±(k) 2 < 4£'|:(k)£ , ±(k). We see 
that for the equality case, we develop a pole for the parameters a±(k), /3ji(k),aq_(k) and /3±(k). In 



an effort to advert this, we consider the condition ae± < y Ae±E ± on the control parameters, where 

e± =el+e b ± . 

Since the injection rates per mode for the kinetic-helicity and magnetic helicity depend on k, the 
global kinetic-helicity injection rate h and the global magnetic-helicity injection rate \ WQ found as, 



at 

dH" 1 



dt 



(el-s u _)J2k = h, (24) 



= (4-^)Er = x- (25) 

p n, 
J S f 



From the above expression of the forces, we find the forces that act in the Elsasser form of the MHD 
equations as f ± = f" ± f . If we chose a purely mechanical forcing of the turbulence (f h = 0) for 
this forcing method, we obtain the relation f + = f = f . For this case e b = 0, which numerically is 
found to be unstable unless no cross-helicity is injected, <r = 0. 

As a note, we see that using the helical decomposition we can redefine the hydrodynamical force we 
used in previous studies [7, 15] (f = 0). From conditions imposed on the energy and kinetic helicity 
injection levels, we take the projections on h± of the force to be: 

f ± (k) = ^. (26) 
N f E ±( k ) 

The condition for the kinetic helicity injection rate is automatically fulfilled for any selection of e±. 



V. CONCLUSIONS AND DISCUSION 



A forcing mechanism which controls the injection level of the three ideal invariants of MHD tur- 
bulence: the total energy, the cross-helicity and the magnetic helicity has been developed and results 
obtained from the implementation of this force into the TURBO solver have been presented. This type 
of force represents a useful tool for the spectral study of turbulence, since it allows to control the level 
to which all three ideal quadratic MHD invariant fluxes relax to. 

As seen from Figure 3, the presence of a non-zero cross-helicity injection level causes a change in 
the two pseudo-energy dissipation rates, which in turn causes a separation of the respective fluxes levels 
once a statistical stationary regime is reached. Looking at the spectra of E + and E~, we observe a 5/3 
scaling. For imbalanced MHD turbulence, we also observe a difference in the pseudo-energy spectra 
levels. However, the scaling exponent tends to remain 5/3. 
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FIG. 3. Spectra (b) and fluxes (a) at the last time point computed in Figure 2, for the Elsasser pseudo-energy. Left 
panels depict balanced turbulence (a — 0; p c = 0) while right panel depict imbalanced turbulence (a = 0.4; 
P c = 0.6). 



By controlling the cross-helicity injection level, this type of force should provide a help in the study 
of solar wind turbulence, which represents a well known case of imbalanced MHD turbulence. Also, 
controlling the magnetic helicity injection level should be of great help in the field of galactic and solar 
dynamo physics. 
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